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REVISITING THE LEAST-SQUARES PROCEDURE FOR 
GRADIENT RECONSTRUCTION ON UNSTRUCTURED 

MESHES* 

Dimitri J. Mavriplis^ 


ABSTRACT 

The accuracy of the least-squares technique for gradient reconstruction on unstruc- 
tured meshes is examined. While least-squares techniques produce accurate results on 
arbitrary isotropic unstructured meshes, serious difficulties exist for highly stretched 
meshes in the presence of surface curvature. In these situations, gradients are typically 
under-estimated by up to an order of magnitude. For vertex-based discretizations 
on triangular and quadrilateral meshes, and cell-centered discretizations on quadri- 
lateral meshes, accuracy can be recovered using an inverse distance weighting in the 
least-squares construction. For cell-centered discretizations on triangles, both the un- 
weighted and weighted least-squares constructions fail to provide suitable gradient 
estimates for highly stretched curved meshes. Good overall flow solution accuracy can 
be retained in spite of poor gradient estimates, due to the presence of flow alignment 
in exactly the same regions where the poor gradient accuracy is observed. However, 
the use of entropy fixes, or the discretization of physical viscous terms based on these 
gradients has the potential for generating large but subtle discretization errors, which 
vanish in regions with no appreciable surface curvature. 

1 INTRODUCTION 

Current-day unstructured mesh aerodynamic production codes rely almost exclusively on for- 
mally second-order accurate discretizations. The two main approaches for achieving second- 
order accuracy involve centrally-differenced convective terms with added artificial dissipation 
[11, 15] and projection-evolution schemes using linearly extrapolated values based on gra- 
dient reconstruction [4, 1]. Other approaches include fluctuation splitting schemes [5], and 
streamwise upwind Petrov-Galerkin schemes [9], although these approaches have not seen 
widespread use in computational aerodynamics. Although artificial dissipation schemes and 
projection evolution schemes have different origins, the final discretizations are closely re- 
lated. Consider the evaluation of a flux at a control volume interface, as depicted in Figure 
1. The projection evolution scheme requires the solution of an approximate Rieman-solver 
at the interface. For example, the often-used Roe Rieman-solver can be written as [19]: 


Fik = F{Ul, Ur) = Fl{u) + Fr{u) + T|A|r ^{ul - Ur) (1) 
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No. NCC-1-02043. 
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Figure 1: Illustration of flux evaluation at control 
volume interface 

where ul and ur represent the value of the flow variables at the left and right sides of the 
control volume interface. For a first-order scheme, these are simply taken as the values at 
the vertices corresponding to the control volume on either side of the face: 

Ul = Ui ( 2 ) 

Ur = Uk (3) 

To obtain second-order accuracy, the left and right states must be obtained by extrapolating 
the control volume values based on a reconstructed gradient. Thus, the second-order accurate 
scheme is obtained using: 


Ul = Ui + \/Ui.fif (4) 

ur = Uk + Vuk-fkf (5) 

where r*/ represents the position vector drawn from vertex i to the center point of the 
control volume interface. This formulation requires the evaluation of the gradients Vu at the 
mesh vertices. These gradients may be evaluated using a Green-Gauss contour integration 
around the vertex-based control volumes, or by taking a least-squares approximation to the 
gradient at each vertex by constructing a tangent plane which best fits the surrounding 
neighboring data in some (weighted) least-squares sense [3, 8]. The Green-Gauss and least- 
squares constructions are outlined in the Appendix. The least-squares construction may 
include weights on the error terms, leading to different gradient approximations for non-linear 
functions. In all cases, the least-squares constructions represent a linear function exactly for 
vertex and cell-centered discretizations on arbitrary mesh types, while the Green-Gauss 
construction represents a linear function exactly only for a vertex-based discretization on 
simplicial elements (triangles or tetrahedra). Various construction techniques for the least- 
squares gradients have been proposed and discussed in the literature. In this work, a Gramm- 
Schmidt construction [3], and a QR decomposition method [2, 8] have been implemented and 
tested. However, very little differences in the computed gradients has been observed between 
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these two construction techniques, while much more important differences due to the choice 
of weights has been found, as will be shown in the paper. 

Artificial dissipation schemes employ a central difference for the convective terms, and 
augment these quantities by a dissipative term which is required for stability. The flux at 
an interface for a first-order accurate artificial dissipation scheme can be written as: 


Fik = Fi{u) + Fk{u) + a{ui - Uk) ( 6 ) 

where a may be a scalar (scalar artificial dissipation) or a matrix (matrix artificial dissipa- 
tion). In the case where o; is a matrix, a natural choice for a, by analogy with equation (1) 
is: 


a = KT\h\T-^ (7) 

where k is a constant to be determined empirically. If k is taken as unity, then the first-order 
accurate matrix dissipation scheme becomes identical to the first-order accurate projection 
evolution scheme. On structured meshes, second-order accurate artificial dissipation schemes 
are obtained by replacing the first difference in equation (6) by a third difference [12]. On 
unstructured meshes, a second-order accurate artificial dissipation flux can be constructed 
as: 


Fik = Fi{u) F Fk{u) + a (Li{u) - Lk{uj) 
where Li(u) represents an undivided Laplacian operator, taken as: 


( 8 ) 


neighbors 

Li{u)= Y. (uk-Ui) (9) 

k=l 

resulting in an artificial dissipation term which is of the same order as a third difference. 
Thus, the second-order accurate matrix dissipation scheme can be obtained by replacing 
the difference of reconstructed states in the projection evolution scheme by a difference of 
undivided Laplacian operators. Although these quantities are of the same order, they are 
not directly proportional to each other, and therefore the parameter k cannot be taken as 
unity in this case, but must be determined empirically. There are also discrepancies be- 
tween the centrally differenced convective fluxes in both schemes, since these are evaluated 
at reconstructed states in the upwind scheme, rather than at vertex values as in the artifi- 
cial dissipation scheme. However, numerical experiments reveal that these differences have 
virtually no effect on solution accuracy in the subsonic and transonic regimes. 

The T matrices on the right hand side of equation (1) represent the eigenvectors associated 
with the linearization of the equations of inviscid compressible flow normal to the control 
volume face ik, while the |A| matrix is a diagonal matrix containing the absolute values of 
the five eigenvalues associated with these equations. Of these five eigenvalues, three are 
repeated, leaving three distinct eigenvalues which are proportional to: u, u+c, u-c, where 
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u is the velocity normal to the control volume face, and c is the speed of sound. When 
one of these eigenvalues vanishes, the dissipation for that component at that location also 
vanishes, which may lead to numerical instabilities. For this reason, it is common to limit 
the eigenvalues to a minimum fraction of the maximum eigenvalue, such as: 


u=sign{u) * max{\u\,S{\u\ + cj) (10) 

u + c=sign{u + c) >i= max{\u + c|, (5(|n| + c)) (11) 

u — c=sign{u — c) * max{\u — c|, (5(|n| + c)) (12) 

where |n| +c is the maximum eigenvalue, and (5 is a parameter to be chosen empirically which 
varies between 0 and 1. When S is taken as 0, no eigenvalue limiting is applied. When S is 
taken as 1, the |A| matrix reverts to a scaled identity matrix, since all eigenvalues are now 
taken as |n| + c, and the triple matrix product r|A|T“^ reduces to a scalar quantity. For 
the artificial dissipation discretization, this constitutes the definition of the scalar artificial 
dissipation, i.e. 


a = K* max eigenvalue 


(13) 


which can be computationally cheaper than requiring the evaluation of the full matrices. 
Small values of S of the order of 0.1 are common in many production codes, and this process 
is often referred to as an entropy fix. 

For fiows with strong gradients, most notably in the vicinity of shock waves, the above 
second-order accurate formulations may lead to instabilities, and additional dissipative mech- 
anisms are required. In the upwind scheme, these take the form of limiters applied to the 
computed gradients [4], while in the artificial dissipation schemes, the differences of undi- 
vided Laplacian operators is replaced by a blend of first differences and undivided Laplacian 
operators [11, 14]. In both cases, accuracy is reduced from second to first order locally 
in regions where this additional dissipation is required. For transonic fiows with shocks of 
moderate strength, the use of limiters or additional dissipation is generally not required. For 
the purposes of the current study, we will confine ourselves to cases where no limiting or 
additional dissipation is employed. 


2 MOTIVATION 

The motivation for the current study comes from the observed behavior of various unstruc- 
tured mesh discretizations for a Reynolds-averaged Navier-Stokes solver on problems of aero- 
dynamic interest. The viscous transonic flow over a wing-body configuration has been com- 
puted using the vertex-based unstructured mesh flow solver NSU3D [17] using a matrix 
artificial dissipation discretization and an upwind scheme based on the unweighted least- 
squares gradient construction technique, and the sensitivity of the solution to the different 
discretizations has been investigated [16]. The particular test case is taken from the 1st 
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AIAA Drag prediction workshop [13]. Figure 2 illustrates the mesh and a sample solution 
for a Mach number of 0.75, a Reynolds number of 3 million and = 0.6. The mesh contains 
a total of 1.65 million vertices, and uses prismatic elements in the boundary layer regions 
and tetrahedral elements elsewhere. The solution is shown as computed surface pressure 
contours, illustrating the shock wave on the upper surface of the wing. 



Figure 2: Computed surface pressure contours 
for transonic flow over wing-body configura- 
tion. Mach=0.75, Cl = 0.6, Reynolds number 
= 3,000,000 

Figure 3 depicts the computed lift and drag values using both the artificial dissipation 
discretization, with an entropy fix value of (5 = 0.1 and the least-squares upwind discretization 
scheme using the entropy fix value of (5 = 0.0 as well as (5 = 0.1 

The lift values produced by the upwind scheme with no entropy fix are slightly lower 
than those given by the matrix dissipation case. Comparing these differences with the 
increased lift values reported by the matrix dissipation on a finer grid of 13 million points, 
one can conclude that the least-squares based discretization is slightly more diffusive than 
the matrix dissipation discretization. Because the nominal value of the k coefficient in the 
matrix dissipation scheme has been determined empirically, it is conceivable that a simple 
rescaling of the dissipation terms could be used to improve the accuracy in the upwind scheme 
as well. However, there are significant differences between these two discretizations which 
extend beyond the simple scaling of the final terms. When the entropy fix parameter for the 
artificial dissipation scheme is increased from 6 = 0.0 to (5 = 0.1, which is the level used in the 
baseline matrix dissipation settings, the results of the upwind scheme are now much different 
than in either baseline cases. The lift is reduced by over 20% and the drag values in the 
polar plot are substantially overpredicted. In essence, the accuracy of the upwind scheme 
has been completely compromised by this small value of the entropy fix. For the matrix 
dissipation scheme, previous studies [16] have shown the solution accuracy to be insensitive 
to small values ((5 = 0.1 to 0.2) of the entropy fix parameter, while the scalar dissipation 
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scheme ((5 = 1.0) achieves reasonable accuracy in lift with slight drag overprediction (si 25 
counts) for the case shown above [16]. The unexpected sensitivity of the upwind scheme to 
small values of the entropy fix prompted the current investigation. 




Figure 3: Variations in computed lift and drag coefficients as a function of discretization 
for transonic flow over wing-body configuration. di denotes the entropy fix parameter 5. 
ALPHA denotes the angle of attack for tha aircraft geometry. 

To better study this problem, we resort to a simpler two-dimensional example. Figure 
4 illustrates the mesh for computations of viscous transonic fiow about an RAE2822 airfoil 
at Mach=0.73, Reynolds=6.5 million, and an incidence of 2.31 degrees. The mesh contains 
a total of 16167 vertices, with the distance of the first point normal to the wall being 2.e- 
06 chords. Although the figure shows a fully triangular mesh, quadrilateral elements are 
employed by the solver in the boundary layer regions by removing the diagonal associated 
with pairs of stretched triangles. A fiow solution behavior similar to that discussed for the 
three-dimensional example is observed in Figure 5. The solutions using the matrix artificial 
dissipation scheme and the least-squares based upwind scheme with a vanishing entropy fix 
agree closely, while the upwind scheme accuracy degrades severely when the value (5 = 0.1 
for the entropy fix is used. 
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Figure 4: Two-dimensional unstructured 

mesh for computation of transonic flow over 
RAE2822 airfoil 



Figure 5: Variations in computed surface pres- 
sure coefficients as a function of discretization 
for two-dimensional transonic airfoil problem. 
Mach=0.73, Incidence=2.31 degrees, Re = 6.5 
million 
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3 INVESTIGATION 


To investigate the cause of this accuracy degradation, we study the accuracy of the various 
gradient construction techniques. Clearly, all the employed gradient construction techniques 
produce the exact result for a linear function. Thus, we must investigate the accuracy of 
these constructions for non-linear functions, and preferably functions which are representa- 
tive of the types of gradients found in real flow simulations. We seek an analytic function for 
which an exact value of the gradient is available, against which the discrete gradient values 
can be compared. This is achieved using the distance function, which represents the distance 
from any given point in the plane to the nearest point on the airfoil surface. Contours of 
the distance function are shown in Figure 6. This represents a convenient choice, since the 
distance function has similar characteristics to boundary layer velocity gradients, (i.e. ex- 
hibits strong normal gradients, and vanishing streamwise gradients) and is readily available, 
since it is required in the formulation of the Spalart-Allmaras turbulence model [20]. By 
deflnition, the gradient of the distance function in the direction normal to the surface is 
unity: 


VT» . n = 1.0 (14) 

and the streamwise gradient vanishes. Thus, the norm of the derivative is given by: 

doV (odV _ ^ 

Figure 7 illustrates the contours of the percent error between the unweighted least-squares 
computed gradient of the distance function and the exact value. For regions away from 
the leading and trailing edges, the maximum error is no more than 0.5%, illustrating the 
accuracy of this construction for nearly linear functions. Similar results are observed with 
the Green-Gauss gradient construction. 

A more non-linear function is constructed by considering a quadratic variation in D: 

F={l + aDf (16) 

At any point, the norm of the gradient in F is given by: 


(15) 


or 


||VF|| = 2a{l + aD) 




||VF|| = 2a{l + aD) 


(17) 


(18) 


following equation (15). For the purposes of this study, the value a = 200 was used ex- 
clusively. This choice was motivated by the desire to minimize roundoff error problems for 
small values of D. 
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Figure 6: Contours of distance function on 
grid of Figure 4 


Figure 7: Contours of percent error be- 

tween computed and exact distance func- 
tion gradient. Errors are below 0.5% ex- 
cept near the leading and trailing edges. 


In Figure 8, the ratio of the computed to the exact value of || VF|| is plotted at the station 
x=0.3 on the airfoil, using the Green-Gauss gradient construction, and the unweighted least- 
squares gradient construction. Additionally, the values obtained from a weighted least- 
squares gradient construction (using inverse distance weighting) are shown, as well as the 
value of dF/dn obtained by finite difference along the normal grid line in the boundary 
layer region. The Green-Gauss and the finite-difference approach produce very accurate 
estimates of the gradient in all regions of the domain. However, the unweighted least-squares 
construction is seen to grossly underpredict the gradient near the wall, and throughout a 
large inner portion of the boundary layer region. When inverse distance weighting is used 
in the least-squares approach, accuracy similar to that achieved by the other methods is 
recovered. Figure 9 shows the value of the gradient at the first grid point away from the 
airfoil surface (in the normal direction), plotted along the entire upper and lower surfaces 
of the airfoil. Once again, the gradient values are well predicted by all methods except the 
unweighted least-squares construction, which shows severe under-prediction and considerable 
scatter aft of the mid-chord location. 
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Figure 8: Ratio of computed to exact gra- 
dient of function F using various methods 
along normal station at x=0.3 location on 
RAE2822 airfoil grid using quadrilateral 
elements in boundary layer region 


Figure 9: Ratio of computed to exact gradi- 
ent of function F at first interior grid point 
closest to wall along RAE2822 airfoil 



Figure 10: Rlustration of unstructured 
mesh employed for computation of vis- 
cous flow over flat plate with rounded 
and tapered leading edge 



Figure 11: Ratio of computed to exact gradi- 
ent of function E at first interior grid point 
closest to wall along upstream portion of flat 
plate geometry 
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To obtain insights into this behavior, we perform the same experiment on a flat plate 
geometry. The grid for this case is shown in Figure 10. The geometry consists of a flat plate 
with a rounded and tapered leading edge. Figure 11 illustrates the comparison between 
analytical and computed gradients at the first point off the wall as a function of streamwise 
coordinate. In this case, the unweighted least-squares gradients compare poorly with the 
exact values near the leading edge, but compare favorably over the downstream region of 
the plate. In fact, the sudden increase in accuracy for these gradients occurs precisely at the 
location where the plate surface becomes horizontal, or more importantly, where the surface 
curvature vanishes. 

This provides an indication as to the mechanism subverting the accuracy of the un- 
weighted least-squares gradient construction. This is illustrated in Figure 12 using the sim- 
plest possible configuration, i.e. a highly stretched quadrilateral mesh in the presence of 
surface curvature. Without loss of generality, we assume the surface normal at the station 
under consideration to be in the y-coordinate direction, and plot the topology using an 
expanded scale in the normal direction. Due to the surface curvature, the upstream and 
downstream neighbors are not aligned with the center point in the y-coordinate. While 
this y-direction variation is indeed very small (order l.e-04 chords on the RAF 2822 mesh 
near the mid-chord location), it is nonetheless much larger than the small normal spacing 
used in the inner portion of the boundary layer region (2.e-06 chords). Therefore, for the 
unweighted least-squares gradient, these points exert a large influence on the determination 
of the normal gradient, in spite of the fact that they are much more distant from the point 
under consideration than the two neighboring values in the upper and lower y-direction. 
This is an unavoidable consequence of the use of an unweighted procedure, which treats all 
(neighboring) stencil points equally. Using the inverse distance weighting in the least-squares 
construction deemphasizes these distant upstream and downstream points, thus resulting in 
much more accurate gradients in such situations. Referring to Figure 12, the accuracy of 
the unweighted least-squares gradient is seen to break down when the normal grid spacing 
h becomes comparable to the distance H. Writing H as a function of the angle a yields 
the necessary condition for avoiding accuracy breakdown of the unweighted least-squares 
gradient: 


h > R{1 — cosa) 


(19) 


or 



( 20 ) 


2 

where s denotes the streamwise grid spacing, and the approximations cosa ~ 1 + ^ + ... 
and ~ have been used. Considering a circle of unit radius, discretized with O(IO^) 
circumferential points, equation (20) predicts poor gradient accuracy for normal grid spacing 
below 0(10"^). This rough estimate correlates well with the behavior observed in Figure 8 
for the airfoil case. 

The failure of the unweighted least-squares gradient construction is perhaps surprising be- 
cause this method possesses many of the often sought-after properties for numerical schemes: 


• It represents a linear function exactly, for arbitrary grid topologies. 
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It has been shown to produce superior gradient estimates for highly irregular (but 
isotropic) meshes. 


• It performs well for cases with no surface curvature, such as flat plate boundary layer 
cases, where most numerical investigations of viscous flow solvers are initiated. 

• It has often been found to be more robust for viscous flows. 

Although the combination of high mesh stretching with surface curvature may be con- 
sidered pathological situations in some disciplines, such topologies are common-place for 
aerodynamic simulations and better gradient estimates are desirable. 



Figure 12: Illustration of stencil for least- 

squares gradient calculation on quadrilateral 
mesh in the presence of stretching and surface 
curvature 


4 EFFECTS WITH ALTERNATE DISCRETIZATIONS 

The above discussion was confined to vertex-based schemes operating on prismatic (3D) or 
quadrilateral (2D) element meshes in the boundary layer region. In this section we examine 
the suitability of the various gradient construction methods for vertex-based discretizations 
on triangular boundary layer meshes, and for cell-centered discretizations using fully trian- 
gular or mixed triangular and quadrilateral meshes. Figure 13 illustrates the topology of the 
least-squares stencil for a vertex-based triangular mesh, and Figure 14 depicts the estimates 
of the gradient of the function F produced by the various methods. In this case, the stencil is 
augmented by two additional points joined by the triangle diagonals, which are at upstream 
and downstream locations from the point under consideration. These additional points are 
similar in character to the upstream and downstream points obtained from the quadrilateral 
stencil, and thus both the unweighted and weighted least-squares methods retain similar per- 
formance on triangular meshes as on quadrilateral meshes. Additionally, the Green-Gauss 
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approach is seen to yield similar results on triangular meshes as in the previous case on 
quadrilateral meshes. Note that the Green-Gauss construction is exact for linear functions 
on triangular meshes only, although this does not appear to have any appreciable effect on 
the accuracy of the results shown in Figure 14. 



Extra Stencil Point 


Figure 13: Illustration of stencil for least- 
squares gradient calculation on triangular 
mesh in the presence of stretching and sur- 
face curvature 



Figure 14: Ratio of computed to exact gra- 
dient of function F using various methods 
along normal station at x=0.3 location on 
fully triangular RAE2822 airfoil grid 


For cell-centered schemes operating on quadrilateral meshes, the stencil is topologically 
similar to that of a vertex scheme operating on a quadrilateral mesh, as shown in Figure 15. 
Thus similar behavior for the various gradient construction methods can be expected. For 
a cell-centered method operating on a fully triangular mesh, the stencil topology is shown 
in Figure 16. In all cases, a stencil with only three neighbors is obtained, and none of these 
neighbors are located close (within the order of a cell normal height) to the cell center under 
consideration. Hence, it is not surprising that the unweighted least-squares gradient exhibits 
poor accuracy in the boundary layer region, similarly to the vertex discretization cases. 
However, inverse distance weighted least-squares construction also exhibits poor accuracy 
in these cases, as shown in Figures 17 and 18, since there are no close points to provide 
accurate normal derivative information. For cell-centered discretizations, the Green-Gauss 
gradient construction generally will not produce the exact value for a linear function. This 
is only achieved if the segments joining neighboring cell centroids bisect the mesh edges, 
which is generally not achieved [4]. For the function F, the gradient values are either over- 
predicted or under predicted by roughly 10 % depending on the orientation of the triangle 
diagonal, as seen by the oscillatory behavior in the plot of Figure 18. (Only one branch of 
these two triangle types is plotted in Figure 17.) The average of these two triangle estimates 
closely approximates the exact gradient value, which is equivalent to performing the integral 
around the quadrilateral formed by the union of the two constituent triangles. In spite of this 
shortcoming on triangular meshes, the Green-Gauss gradients are seen to provide superior 
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estimates of the gradients of F in the boundary layer regions to the least-squares methods. 
This illustrates the danger of relying on simple properties such as exact representation of 
linear functions for accuracy certification. 




Figure 15: Illustration of stencil for 

least-squares gradient calculation for cell- 
centered discretization on quadrilaterals 


Figure 16: Illustration of stencil for 

least-squares gradient calculation on cell- 
centered discretization on triangular mesh 
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Figure 17: Ratio of computed to exact gra- 
dient of function F along normal station at 
x=0.3 location for cell-centered discretiza- 
tion on fully triangular RAE2822 airfoil 
grid 


Figure 18: Ratio of computed to exact gra- 
dient of function F at second layer of trian- 
gles off wall for cell-centered discretization 
on fully triangular RAE2822 airfoil grid 
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Alternate techniques for constructing gradients on cell centered simplicial discretizations 
have been developed [6, 10, 7, 18]. In one of these approaches [7, 18], vertex values are 
obtained by averaging surrounding cell-centroidal values, often using a weighting factor, and 
the cell based gradient is then computed using a Green-Gauss contour integral using the con- 
structed vertex values. The technique described in [6, 10] extrapolates gradients computed 
on neighboring triangles using a Green-Gauss contour integration. The performance of these 
strategies for highly stretched meshes in the presence of curvature has not been studied in 
this work, but warrants further investigation. 

5 EFFECT ON SOLUTION ACCURACY 

While we have pointed out inadequacies in the unweighted least-squares gradient formulation, 
the fact remains that in the absence of any entropy fix, upwind discretizations based on this 
approach achieve good overall accuracy, as evidenced by the results in Figures 3 and 5. It 
may seem perplexing how one can obtain a viable solution with such poor gradient estimates 
in the inner boundary layer region, and why this solution is so sensitive to small values of the 
entropy fix parameter. The answer lies in the alignment of the grid with the flow direction 
in the boundary layer region. The use of a highly stretched mesh aligned with the wall 
direction in boundary layer regions, which is commonplace for high-Reynolds number fow 
simulations, results in near vanishing flow velocity normal to the control volume interfaces 
in this direction. This is shown in Figure 19, where the computed normal and streamwise 
flow velocities are plotted along the normal station at x=0.3 for the transonic airfoil flow 
solution using matrix dissipation depicted in Figure 5. This plot indicates that the normal 
velocity is two to three orders of magnitude smaller than the streamwise velocity throughout 
the inner portion of the boundary layer. In Figure 20, the normal and streamwise convective 
eigenvalues are plotted at the same station. The convective eigenvalue is defined as the 
integrated velocity flux through the control volume face. Due to the high cell aspect ratios, 
the normal control volume face is much larger than the streamwise face, and the normal 
eigenvalue becomes substantially larger than the streamwise value in the inner boundary 
layer region. Thus, in spite of decoupling through flow alignment, the high grid cell aspect- 
ratio ensures strong coupling between neighboring normal points in the boundary layer, even 
for the convective modes. However, the overall dissipation terms are formed by the product 
of the eigenvalue with the jump in left and right flow variables across the control volume face. 
Assuming (as a worst case scenario) a first order variation in the flow variables, the normal 
(streamwise) dissipation terms scale as the product of the normal (streamwise) eigenvalue 
and the normal (streamwise) grid spacing, with this latter quantity being of the same order 
as the streamwise (normal) control volume face. Thus the overall scaling of the dissipation 
terms is closely approximated by Figure 19, which implies much lower diffusion in the normal 
direction as compared to the streamwise direction. The application of an entropy fix places 
a lower limit on the velocity values used in scaling the dissipation terms, which from Figure 
19 can be seen to have a large effect on the normal velocity values. 

It is interesting to note that, although the dissipation terms associated with the normal 
velocity are small, the two acoustic wave eigenvalues associated with u+c and u-c are not 
affected by flow alignment, and yet good accuracy is retained despite the use of inaccurate 
gradients for the dissipative terms associated with the acoustic waves. On the other hand, the 
use of more accurate gradient estimates resolves the loss of accuracy for small values of the 
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entropy fix parameter. In Figure 21, the transonic airfoil fiow case has been recomputed using 
the upwind scheme with an inverse-distance weighted least-squares gradient construction, 
using a vanishing entropy fix, as well as an entropy fix parameter value of (5 = 0.1. The 
computed surface pressures in both cases compare well with each other and agree closely 
with those produced by the matrix artificial dissipation scheme, illustrating the superior 
characteristics of the weighted versus the unweighted least-squares construction. 



Figure 19: Magnitudes of normal and 

streamwise velocities along normal station 
at x=0.3 for transonic airfoil flow solution 



Figure 20: Magnitudes of normal and 

streamwise convective eigenvalues along 
normal station at x=0.3 for transonic air- 
foil flow solution 
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Figure 21: Variations in computed surface 

pressure coefficients for weighted least-squares 
gradient construction as a function of dis- 
cretization for two-dimensional transonic air- 
foil problem. Mach=0.73, Incidence=2.31 de- 
grees, Re = 6.5 million 


6 IMPLICATIONS 

In the above discussion, we have demonstrated how and why the unweighted least-squares 
gradient construction severely under-estimates normal gradients for highly stretched meshes 
in the presence of surface curvature. Furthermore, it has been shown that this behavior 
can be expected for vertex based discretizations operating on triangular and quadrilateral 
meshes (tetrahedral and prismatic meshes in 3D) and for cell centered discretizations on 
either types of meshes as well. The use of inverse distance weighting in the least-squares 
construction can be used to recover good accuracy in these situations for vertex and cell 
centered discretizations on quadrilateral meshes, and for vertex discretizations on triangular 
meshes. However, this technique is not effective for cell centered discretizations on triangular 
meshes. The Green-Gauss construction technique produces adequate gradient estimates in 
all cases, even for cell centered discretizations where it may not represent linear functions 
exactly. 

On the other hand, the failure of the un-weighted least-squares gradient is mitigated by 
the phenomenon of flow alignment in precisely the same locations, thus enabling adequate 
overall accuracy to be achieved. Therefore, the least-squares gradient construction can be 
used competitively for producing accurate solutions, but the user must be aware of the 
limitations of this approach: notably that no entropy fix be used, and that the mesh be 
well aligned with the viscous surfaces, and thus with the boundary layer flow direction. 
Alternately, mesh cell aspect ratio constraints based on equation (20) may be enforced. 
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Discretization of the physical viscous terms in the Navier-Stokes equations is often 
achieved in a two pass approach where the flow gradients are evaluated in the first pass, 
and then used in the second pass to build up these terms. Clearly, the use of unweighted 
least-squares gradients in the construction of these viscous terms has the potential to gen- 
erate large discretization errors overall. However, the deficiency of this approach may be 
extremely subtle in that it will not at all manifest itself for flat plate boundary layer calcu- 
lations, where most viscous flow solvers are initially validated, but only in the presence of 
bodies with non-negligible surface curvature. 

Finally, a more prudent strategy would appear to be one which employs inverse distance 
weighted least-squares gradients or even Green-Gauss gradients in the discretization of con- 
vective and viscous terms. However, these approaches have proven to be substantially less 
robust than upwind schemes based on unweighted least-squares gradients, and often require 
gradient limiting to achieve stable solutions, which in turn, may have an adverse impact on 
accuracy. While it was initially argued that this was due to superior approximation proper- 
ties of the unweighted least-squares approach especially for irregular meshes, it should now 
be evident that the main reason for the robustness of this approach can be attributed to the 
use of under-predicted gradients, effectively using limited gradients which correspond to a 
first order scheme in the inner part of the boundary layer. 

An alternative to the inadequacies of the unweighted approach, and the poor robustness of 
the weighted approach, is to resort to different gradient formulations, such as those described 
in [18, 7, 6, 10], or to employ one-dimensional reconstruction and limiting in the normal 
direction in the boundary layer, analogous to structured mesh techniques, although this 
incurs obvious data-structure drawbacks for an unstructured mesh approach. Future work 
will investigate the accuracy and robustness of various such discretizations for both vertex- 
based and cell centered approaches. 

7 APPENDIX 

The Green-Gauss formulation constructs gradients by integrating around the boundary of a 
closed control-volume. From Green’s theorem, the average gradients over a control volume 
can be written as: 




{Uy)i 


1 


Vok 

1 


Vok 


UxdV 

UydV 



( 21 ) 

( 22 ) 


where Vok represents the area of the two-dimensional control volume. These contour inte- 
grals are then discretized as: 


' ~ Vok £ 2 

(23) 

Vok ti 2 

(24) 


where the x and y subscripts denote differentiation, and i and k identify the associated 
control volume. For vertex schemes, the median dual control volumes are employed, as 
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depicted in Figure 1. In this case, i and k refer to the vertices on either side of the control 
volume face, and Axik and Ai/ik denote the increments of x and y along the control volume 
face. For cell-centered schemes, the grid cells themselves form the control volumes. In this 
case, i and k refer to the cells on either side of a mesh edge, and Axik and Ayik denote 
the increments of x and y along the mesh edge. The Green-Gauss formulation is exact 
for vertex discretizations of linear functions only on triangular elements. For cell-centered 
discretizations, this formulation is generally not exact for linear functions on quadrilaterals 
or triangles. In the special case where the segments joining neighboring triangle cell centers 
exactly bisect the shared mesh edge, the formulation becomes exact for linear functions on 
triangles. 

The least-squares gradient construction is a technique which is unrelated to the mesh 
topology. This construction relies on a stencil which identifies relevant neighboring points 
for use in the gradient estimation. Although this stencil can be chosen arbitrarily, the 
most obvious construction for mesh-based data is to chose the stencil of nearest neighboring 
values. For a vertex-based discretization, the stencil is thus formed by the set of mesh edges 
incident on the considered vertex i. For a cell-centered discretization, the stencil is formed 
by the edges joining neighboring cell centroids, which corresponds to the dual graph of the 
mesh, as shown in Figures 15 and 16. The least-squares gradient construction is obtained 
by solving for the values of the gradients which minimize the sum of the squares of the 
differences between neighboring values k = 1,N and values extrapolated from the point i 
under consideration to the neighboring locations. The objective to be minimized is given as 

N 

E (25) 

k=l 


The error is given as 


^ik ( dUiP; + (^Ux)i-dXpp; + i'^y)i-dyik) 


(26) 


where duik represents the difference Uk — Ui, with similar expressions for dxik and dyik, and 
Wik is a weighting factor. Dropping the i subscripts on the gradients for clarity, a system 
of two equations for the two gradients and Uy is obtained by solving the minimization 
problem i.e. 




ik 


dux 


= 0 


(27) 


and 


^lE. 


ik 


du„ 


= 0 


(28) 


Using straight-forward algebra, the following set of equations is obtained 


^i'^X billy di 

billx 


(29) 

(30) 
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where 


0,1 


(31) 

h = 

T,k=i wf^dxikdyik 

(32) 

Ci = 

Eti^Ikdyfk 

(33) 

di 

T,k=i wj^duikdxik 

(34) 

= 

T,k=i wf^duikdyik 

(35) 


In practice, all the above terms can be precomputed and stored, since these are only a 
function of the grid metrics. The above system of equations for the gradients is then easily 
solved using Cramer’s rule. Note that the determinant of this system is given by: 

DET = ac — (36) 

For the unweighted case {wik = 1), the determinant corresponds to a difference in quantities 

of the order 0{dx'^), which may lead to ill-conditioned systems. This may be the motivation 

for investigations into alternate solution techniques for the least-squares construction, such 

as the QR factorization method advocated in [2, 8]. For the flows computed in this work, 

very little difference in the calculated gradients was observed between these methods. Note 

that when inverse distance weighting is used (wik = , =), the determinant scales as 

V'^^ik+'^yik 

0(1), and the system is much better conditioned. 
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